Monte Carlo Study of an Inhomogeneous Blume-Capel Model 
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Systems of particles in a confining potential exhibit a spatially dependent density which fundamen- 
tally alters the nature of phase transitions that occur. A specific instance of this situation, which is 
being extensively explored currently, concerns the properties of ultra-cold, optically trapped atoms. 
Of interest is how the superfluid-insulator transition is modified by the inhomogeneity, and, indeed, 
the extent to which a sharp transition survives at all. This paper explores a classical analog of these 
systems, the Blume-Capel model with a spatially varying single ion anisotropy and/or temperature 
gradient. We present results both for the nature of the critical properties and for the validity of the 
"local density approximation" which is often used to model the inhomogeneous case. We compare 
situations when the underlying uniform transition is first and second order. 

PACS numbers: 71.15.Mb, 37.10.Jk, 64.60.De, 64.60. fd 



Introduction 

The realization of superfluid and Mott insulator tran- 
sitions in optically trapped atomsi*^!^ has led to an ex- 
amination of the nature of phase transitions in the pres- 
ence of a spatially varying potential. For example, it 
was found that when a confining potential is added to 
the Bose-Hubbard Hamiltonian, the variation of density 
across the sample results in a coexistence of superfluid 
and Mott insulator regions^i^. As a consequence, crit- 
ical phenomena which occur in the uniform case, when 
the entire system collectively makes a transition from one 
phase to another, are smeared. The density no longer ex- 
hibits a singularity as a function of chemical potential, as 
occurs in the translationally invariant caso^. Measures 
of "local quantum criticality" can be defined to help draw 
out residual signals of the global phase transition^. 

These conclusions have been drawn from direct ex- 
amination of the inhomogeneous model, but have also 
been inferred from studies of the translationally invari- 
ant model combined with the "local density approxima- 
tion" (LDA)^. Specifically, the LDA assumes that the 
properties of the confined system at a particular spatial 
location are identical to those of the unconfined system 
with a uniform potential taking the same value as the lo- 
cal potential at that location. Various checks have been 
made, for example by comparing the LDA results using 
quantum monte carlo (QMC) of a collection of uniform 
systems, with QMC simulations of a lattice with a real 
trap^i^. 

This LDA approximation is of course in direct analogy 
with that commonly used in density functional theoryiS, 
where the exact exchange-correlation potential present 
at a particular position r, in a system where the electron 
density varies spatially, is replaced by the exchange cor- 
relation energy of the uniform electron gas at the same 
constant density as that present at r. It is known that 
this approximation yields very good results in a number 
of contexts, especially when the electron-electron interac- 
tions are of weak to intermediate strength. On the other 
hand, when the coupling is stronger, and phenomena like 



magnetism and Mott transitions occur, the LDA is less 
accuratoii. 

In this paper, we examine the nature of phase transi- 
tions in spatially inhomogeneous systems, and the valid- 
ity of the LDA, within a more simple classical context. 
Previous work in this area includes studies of Ising tran- 
sitions in systems with a temperature gradient where the 
nature of the interface between ferromagnetic regions ad- 
jacent to the "cold side" (T < Tc) of a sample and para- 
magnetic regions next to the "hot side" (T > Tf.) has 
been exploredi^ii^i^ii^. 



Model and Calculational Approach 

A classical model which can be constructed to have 
a spatially varying density similar to that in optically 
trapped atom systems is the Blume-Capel modeli^iii 
with a site dependent single-ion anisotropy, 
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Here Si is a discrete classical variable which can take on 
three values. Si = 0, ±1. A coupling J is present between 
near-neighbor spins which we choose to be positive (fer- 
romagnetic). We consider a square lattice of linear size 
L. That is, i = (ix,iy) with 1 < ix,iy < L. The value 
Si ~ Q can be thought of as corresponding to a vacancy, 
while 5*1 = ±1 is an Ising spin, a collection of which can 
order ferromagnetically if the ratio of J to temperature 
T is sufficiently large. A is the single-ion anisotropy pa- 
rameter, and controls the density of Si = spins. 

The Blume-Capel model was originally introduced by 
Blume>i^ and Capelii, separately, to study first-order 
magnetic transitions. It was later generalized to the 
Blume-Emery-GrifBths model (BEG)^, which incorpo- 
rates an additional biquadratic interaction K SfS?. 
Since their initial fomulation, the Blume-Capel and BEG 
models have been extensively used to study the phase 
separation of He'^-He* mixtures^^ and various other sys- 
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terns that exhibit tricritical behavior, such as multi- 
component fluidsi^ and semiconductor alloys2S. Re- 
cent works have used the Blume-Capel model to study 
ferromagnetic thin films using an alternating singie-ion 
anisotrop}*2i, and the dynamics of rough surfaces^. 

Our computational method is standard Metropolis 
Monte Carlo. Each spin of the lattice is visited and a 
change from the current spin value to one of the two other 
possibilities is suggested. This change is accepted or re- 
jected with the Metropolis prescription. To ensure equi- 
libration, a large number of sweeps of all the spins in the 
lattice is performed prior to making measurements. Un- 
less otherwise noted, the statistical errors in our results 
are smaller than the symbol size. The lattices studied in 
the manuscript are small enough that it is not necessary 
to emply more powerful cluster algorithms such as those 
developed by Swendsen and Wang2^. 

For uniform systems, an accurate determination of the 
critical point can be obtained from computing the second 
moment of the magnetizatiori^i. 



(2) 



Near the critical temperature, Tc, the following finite size 
scaling expression holds. 



(3) 



Here (3 {v) are the critical exponents governing how the 
magnetization (correlation length) vanishes (diverges) as 
T Tc in the thermodynamic limit. Eq. 3 implies that 
plots of L^f^/" (AP) for different lattice sizes L cross at 
T ~ Tc, providing a method to locate the critical tem- 
perature. 

The physics of the Blume-Capel model with uniform 
Ai = A is well understood. When A — > — oo, vacancies 
{Si = 0) are energetically very unfavorable. The system 
reduces to the Ising model and there is, on a square lat- 
tice, a second order magnetic phase transition at Tc = 
2.269J. We can also deduce the critical coupling at zero 
temperature. The energy of the fully polarized ferromag- 
netic state (all Si = +1) is £:fcrro = (-2 J + A)L'^. The 
energy of the empty state (all S'i — 0) is i?vacuum = 0. 
The ferromagnetic phase is favored up until A > 2J. 
Thus the phase diagram in the (T/J, A/J) plane con- 
sists of a ferromagnetic region at low T/J and low A/J 
bounded by the hnes T/J = 2.269 and A/J = 2. As A 
increases from A = — oo the extra entropy of vacancies 
reduces Tc until the Ising limit boundary bends over to 
contact the T = critical point. 

The phase boundary for uniform J, A has been ob- 
tained by a number of methods, including Monte Carlo 
siinulations2^i^i21, finite-size scalin g^^i^^ , renormaliza- 
tion group methods^i^, and series expansions^. From 
these studies it is known^^ that there is a tricriti- 
cal point along the phase boundary at (T/J, A/J) = 
(0.609(4), 1.965(5)). At low temperatures in the vicin- 
ity of the T = critical point (T/J, A/J) = (0,2) the 
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FIG. 1: (color online) Top: Phase diagram of the Blume Capel 
model at uniform Ai = A. Second (first) order phase tran- 
sitions are indicated by the solid (dashed) lines. The dot- 
ted line is the Ising limit. The four diamonds depict values 
from Ref. [2^. The squares are taken from Ref. l34 . The rest 
of boundary was obtained using our code and the finite size 
scaling analysis of Eq. 3. The arrows denote the trajectories 
used in the simulations of the inhomogeneous system, and 
correspond to Figs. 2,3,4,5 as indicated. See text. Bottom: 
A representative finite-size scaling analysis is shown. Here 
A = 0. The critical temperature Tc is determined by the po- 
sition of the universal crossing of the scaled second moment 
of magnetization for different linear lattice sizes. 



magnetization jumps discontinuously upon leaving the 
ferromagnetic phase. Beyond the tricritical point, the 
transition becomes continuous (second order). The phase 
boundary for this model is shown in Figure [TJtop), where 
the values for Tc were obtained through the analysis of 
Eqs. 2,3 using our code and from Ref. [H andHj. Fig- 
ure [ijbottom) shows a representative finite size scaling 
crossing for A/J = 0. Table U provides the locations of 
Tc for various values of A. 

Having reviewed and reproduced some of the features 
of the transitionally invariant Blume-Capel model, we 
now turn to the subject of this paper, the inhomogeneous 
case. We choose three models of spatial inhomogencity. 
In the first two we introduce a linear variation of either 
the single-ion anisotropy or the temperature, keeping the 
other parameters fixed. 
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TABLE I: Table of the critical temperatures for various values of A/J that were found using our code and finite-size scaling 
teclmique, in comparison with those from Refs.[23 and[29l. The value at A/J — —8, where vacancies are strongly suppressed, 
is close to the Tc/J = 2.269 of tlie two dimensional Ising model, as expected. 



r(i) = To+ ^\ ^° zx, A = const (4) 

These correspond to vertical (A varying) and horizontal 
(T varying) cuts in the phase diagram. In the third case 
we allow both temperature and single ion anisotropy to 
change together, 

A , Ai-Aq . 

A(i) = Ao H ix, 



r(i) = m( Ai + 



(5) 



where m determines the slope of T(i). This more general 
inhomogeneity aUows us to foUow paths in the (T, A) 
plane which are perpendicular to the phase boundary 
in the intermediate coupling regime where the boundary 
curves around from its low T and large negative A limits. 
Typically we will be interested in cases where Aq, Ai,To 
and Ti are chosen such that the lattice is ferromagnetic 
on the left side, ix = 1, with very few vacancies, and then 
becomes paramagnetic for ix = Lx- For simplicity, we 
have chosen a gradient only in one spatial direction x, so 
the iso-contours of the single ion anisotropy are vertical 
lines. In d = 2 ultracold trapped gases, the iso-contours 
are typically circles around the trap center. However, 
we do not expect the results of our study to depend on 
the shape of the boundary between phases, only on the 
existence of the boundary itself. 

We have imposed periodic boundary conditions (pbc) 
in both the x and y directions. Besides reducing finite 
size effects, the use of pbc avoids having edge sites with a 
smaller number of neighbors than in the bulk, a situation 
which would make the connection with the LDA less sim- 
ple. However, there is one slightly tricky issue with the 
pbc. The pbc links in the y direction by construction con- 
nect sites with the same Ai. In the x-direction, the pbc's 
link sites with vastly different single ion anisotropics: Aq 
and Ai. To avoid this problem, the simulations were run 
on lattices with linear size 2Lx + 1 in the x-direction with 
Ai symmetric across the center of the lattice. In effect, 
a second copy of the lattice is connected to the x = \ 
boundary of the first, and the values of Ai increase lin- 
early back up to Ai at which point the pbc connection 



is established. A final point about the geometry is that 
when we explore finite sizes effects we will fix Ly ~ 50 
and increase Lx at constant (Ai — Aq). This is done be- 
cause the X direction is the one along which the gradient 
is established and so increasing Lx allows us to explore 
the limit where the anisotropy gradient becomes weaker 
and weaker. Each of the cases will be used in regions 
where the respective gradient is approximately perpen- 
dicular to the phase boundary, as shown in figure [Ijtop). 
This ensures that the critical region will be localized to 
a small area of the lattice. 

We present results for the "linear structure factor," 
which we define as, 
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S{ix) sums up the spin-spin correlations for all separa- 
tions ly = 1, 2, • • • i with a given ix- The pairs of sites 
in S{ix) therefore all have the same value of Ai. This is 
a convenient (indeed essential) choice in order to make 
meaningful comparisons with the LDA which employs 
lattices of constant A. In this way, S(ix) is the natural 
generalization of the mean square magnetization (Eq. 2) 
used in the translationally invariant case. 

We will also compare the energy for an inhomogeneous 
lattice with that obtained by the LDA. Similar consider- 
ations apply here as with the linear structure factor; we 
would like to compare observables for sets of sites with 
the same value of Ai. However, the energy involves links 
(in the x-direction) which connect sites with different Ai. 
For this reason we will present results for the energy as- 
sociated with bonds only in the y direction. 
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FIG. 2: (color online) The full lines depict the linear struc- 
ture factor S versus A(i) at T/J = 1.25 for several differ- 
ent gradients in the single-ion potential. For all full lines 
Ao/J = —8.00 and Ai/J — 4.00, but as Lx increases the 
gradient 5 A = (Ai — Ao)/Lx softens. The LDA result is the 
dashed curve, and is quantitively correct except in the transi- 
tion region. As expected, the accuracy of the LDA improves 
as the gradient of the inhomogeneity decreases. 



FIG. 3: (color online) S versus A(i) at T/J = 0.56, a tempera- 
ture below the tricritical point. Data are shown for several dif- 
ferent values of the single-ion potential gradient (full curves). 
While the LDA predicts a first order phase transition (dashed 
curve), as is expected at this temperature, the phase transi- 
tions of the inhomogenous systems are less abrupt. As the 
variations in the single-ion potential decrease, or one moves 
far from the transition region, S converges to the results ob- 
tained by the LDA. 



Results 

Figure [5] shows, for fixed T/ J ^ 1.25, the hnear struc- 
ture factor 5 as a function of A/J. More precisely, S{ix) 
is computed at different values of ix for a system with a 
gradient i5A = (Ai - Ao)/Lx with Aq/J = -8.00 and 
Ai/J = 4.00. The value of S{ix) is plotted against the 
corresponding value of A(j^ on the horizontal axis. 
Since the relation between A(j^ and ix is linear (Eq. 4), 
the horizontal axis can equivalently be viewed as labeling 
the spatial position as one sweeps across the inhomoge- 
neous lattice. At the same time, the LDA values are 
obtained by simulating uniform systems at a range of A 
values corresponding to the vertical trajectory marked 
"Fig. 2" in the phase diagram of the uniform system. 
Fig. [TJtop). This trajectory crosses the ferromagnetic to 
paramagnetic phase boundary at A = 1.300(3) in a sec- 
ond order transition. We see that the LDA predicts the 
behavior of iS in a qualitatively correct fashion over the 
entire range of A, and is quantitively accurate except in 
the vicinity of the critical region where the lattice inho- 
mogeneity blurs the transition. This is the same basic 
result as found for optically trapped atom systems^i^. 
However, we are able in this simple classical model to 
compare more precisely the LDA with the inhomoge- 
neous case. In particular. Fig. [2] shows the improved 
accuracy of the LDA as the gradient of the inhomogene- 
ity becomes smaller, something which has not yet been 
done in the quantum case. 

Fig. [3] shows a similar set of data but for T/J = 0.56 
which corresponds to the trajectory labeled "Fig. 3" in 
Fig. 1 and crosses the ferromagnetic-paramagnetic phase 
boundary in a first order transition at A/J — 1.979(2). 
Again, the LDA is qualitatively correct. Comparing 
Figs. [5] and [31 it appears that the LDA has larger quanti- 



tative errors in the vicinity of the transition region in the 
first order case, but that these errors extend less far away 
from the transition region. This result seems reasonable: 
a smoothly varying potential does not exhibit the very 
abrupt discontinuity in the LDA results, but because the 
first order transition region is narrower, the region where 
LDA fails significantly is less wide. It is notable that 
curves of the linear structure factor for different values 
of the gradient cross at roughly a single point. 

Our final two results for the linear structure factor are 
given in Figs. S] and [5l and show cuts across the phase 
boundary in which both temperature and the single ion 
anisotropy are simultaneously evolving. The trajectories 
are labeled Figs. 2] and [SI in Fig. [1] The physics of our 
model does not depend independently on T, Aj, and J, 
but only on the ratios Ai/T and J/T. In Fig.[2|only the 
first of these ratios is changing, while in Figs. [4] and [5] 
both ratios are evolving as we traverse the lattice. Since 
both of these cuts traverse the phase boundary in the 
second order region, the results for the LDA resemble 
those of Fig. [2] This emphasizes that the question of the 
accuracy of the LDA does not appear crucially to depend 
on which parameters in the energy (or the temperature) 
are varying. 

We now turn to a comparison of the energy of the in- 
homogeneous system with that of the LDA. Fig. [Hi shows 
the same cut at constant T/J = 1.25 as in Fig. [51 Re- 
markably, the energy is given very accurately by the LDA 
throughout the inhomogeneous lattice, even through the 
transition region where the linear structure factor dif- 
fered markedly. This result may appear surprising in 
that the first piece of Ey in Eq. 5 is the near neighbor 
spin correlation in the y-direction, which is also one of 
the ingredients of the linear structure factor S. That the 
LDA value for Ey is so accurate suggests that the fail- 
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FIG. 4: (color online) The linear structure factor 5 versus 
r(A(i)) for the inhomogenous case with both a spatially vary- 
ing single-ion potential and temperature gradient. The choice 
m = 0.1458 in Eq. 5 makes this trajectory cross perpendicular 
to the phase boundary. 



1 

<s> 



0.8 



0.6 



0.4 











A/J=-4 












• 50 X 60, ST = 0.0417 






■ 50 X 100, 5T = 0.025 






50 X 300, 5T = 0.0058 






-- LDA 

















1.5 



T/J 



FIG. 5: (color online) The linear structure factor 5 versus 
r(i) at A/J = —4 for several different values of temperature 
gradients. 



ure of the LDA in the transition region is dominated by 
its mis-estimate of the long-range correlations, while the 
short range-ones are correctly identified. Indeed, this re- 
sult might be expected since it is the long-range pieces of 
iS whose behavior is crucial to the occurrence of a second 
order transition. 

Finally, Fig. [7] shows the same cut at constant T/J = 
0.56 as in Fig. [3] Here, when the underlying homoge- 
neous transition is first order, we see that even the energy 
is badly estimated by the LDA. The local energy has a 
universal crossing, corresponding to the transition value 
of A/J, similar to that of the linear structure factor S. 



Conclusions 

The "Local Density Approximation" is a commonly 
employed method to understand the phase transitions 
of ultracold, optically trapped atoms which experience 
a spatially varying confining potential. In this paper we 
have explored the validity of the LDA in the simpler, clas- 
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FIG. 6: (color online) Comparison of the LDA prediction for 
the energy with the energy of an inhomogeneous system. Here 
we have fixed T/J — 1.25 and are changing Ai across the 
lattice. This is the trajectory labeled Fig. [2] in Fig. [T] The 
LDA energy is remarkably accurate even in the transition 
region. 
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FIG. 7: (color online) Comparison of the LDA prediction for 
the energy with the energy of an inhomogeneous system. Here 
we have fixed T/J = 0.56 and are changing A, across the 
lattice. This is the trajectory labeled Fig. [3] in Fig. [T] Near 
the transition the LDA energy differs markedly from that of 
an inhomogeneous lattice. 



sical, Blume-Capel model to which we have applied a gra- 
dient in the temperature and/or the single ion anisotropy. 

Our basic conclusion is that the LDA performs well 
quantitatively in regions that are not close to where the 
state of the system is making a transition between the 
allowed phases, in our case ferromagnetic and paramag- 
netic. That is, the values of the local structure factor and 
energy predicted by the LDA match those of a direct sim- 
ulation of the inhomogeneous system except in the tran- 
sition zone. This is similar to the conclusions drawn in 
the optical lattice case^. However, because our model is 
classical as opposed to quantum mechanical, we can ex- 
plore the validity of the LDA in greater detail, including 
the systematic improvement in the accuracy of the LDA 
with inhomogeneous systems which have smaller gradi- 
ents. 

An especially interesting feature of the Blume-Capel 
model is the presence of a tricritical point on the phase 
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boundary. This allows us to compare the validity of the 
LDA for first and second order transitions in the same 
model. Our conclusion is that a quantity like the linear 
structure factor which samples long range correlations 
is more badly estimated by the LDA in the transition 
region of a first order phase change, but that the width 
of the region over which the LDA is inaccurate is more 
narrow. Overall, the accuracy of the LDA in the two 
cases is not so dramatically different. On the other hand, 
the predictive accuracy of the LDA for the energy, which 



samples just short range correlations, is very different 
for the first and second order situations. In the second 
order case, the LDA energy is quantitativly correct even 
through the transition region, while in the first order case 
the energy is rather badly mis-estimated. 

This work was supported by CNRS (France) PICS 
18796 and ARO Award W911NF0710576 with funds from 
the DARPA OLE Program. We acknowledge useful input 
from T. Hollies. 
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